module equilibrium
   !
   ! Subroutines:
   !
   ! 1) eq_frozen_cp_gamma
   ! 2) eq_effective_gamma
   ! 3) eq_Mach
   ! 4) eq_frozen_mi_k
   ! 5) eq_entropy_and_corners_correction
   !
   ! Last update: 10 Feb 2012

   use portlib
   use Gibbs_1p4

contains

   !****************************************************************************

   ! Subroutine 01

   subroutine eq_frozen_cp_gamma(nx, ny, T, p, cmod, Nr, itt, nitm, itimax, &
         itemax, tol_ej, tol_n, OF, Yi, hi, Rgp, cp, gcp, ej) ! The last three are output

      integer, intent(in) :: nx     ! Number of control volumes in csi direction (real + ghost)
      integer, intent(in) :: ny     ! Number of control volumes in eta direction (real + ghost)
      integer, intent(in) :: cmod   ! Chemical model
      integer, intent(in) :: itimax ! Maximum number of iterations to solve e(j) - dissociation rate for j-reaction
      real(8), intent(in) :: tol_ej ! Tolerance for e(j) evaluation
      integer, intent(in) :: itemax ! Maximum number of iterations to solve n - the number of mols of products
      real(8), intent(in) :: tol_n  ! Tolerance for n evaluation
      real(8), intent(in) :: OF     ! Oxidant/Fuel ratio
      integer, intent(in) :: Nr     ! Number of chemical reactions of the adopted scheme
      integer, intent(in) :: itt    ! Main iteration cycle
      integer, intent(in) :: nitm   ! Maximum number of iterations in main cycle

      real(8), dimension(:), intent(in) :: T     ! Temperature (K)
      real(8), dimension(:), intent(in) :: p     ! Pressure (Pa)
      real(8), dimension(:), intent(out) :: Rgp   ! Mixture of gases constant (J/kg.K)
      real(8), dimension(:,:), intent(out) :: Yi    ! Mass fraction (non-dim.)
      real(8), dimension(:,:), intent(out) :: hi    ! Enthalpy (J/kg.K)
      real(8), dimension(:), intent(out) :: cp    ! Specific heat at const pressure (J/kg.K)
      real(8), dimension(:), intent(out) :: gcp   ! gcp = gamma = Cp/Cv at the center of CV P (non-dim.)
      real(8), dimension(:,:), intent(inout) :: ej    ! chemical reaction rates

      real(8), dimension(Nr) :: e_j  ! auxiliar for chemical reactions

      integer :: nps           ! Counter: south
      integer :: i, j, ii, np  ! Counters


      ! Informations which need to be passed to Gibbs module
      call GIBBS_MACH1D_recebe4 (itt)
      call GIBBS_MACH1D_recebe5 (n)
      call GIBBS_MACH1D_recebe7 (nitm)
      call GIBBS_MACH1D_alocacao_ej

      ! Internal control volumes

      do i = 2, nx-1
         do j = 1, ny

            np = (j-1)*nx + i

            ! Information about the P CV to Gibbs module
            call GIBBS_MACH1D_recebe3 (np)

            call GIBBS_EQUILIBRIO_dados_recebe (1, 0, cmod, itimax, itemax, &
               tol_ej, tol_n, T(np), p(np)/101325.0d0, OF)

            if (itt > 1) call GIBBS_MACH1D_recebe6 (ej)

            call GIBBS_EQUILIBRIO_calculos1

            call GIBBS_MACH1D_envia8 (e_j)

            do ii = 1, Nr
               ej (np,ii) = e_j(ii)
            end do

            call GIBBS_hi_massa_calculo
            call GIBBS_MACH1D_envia4(cp(np), gcp(np), Rgp(np), Yi(np,:), &
               hi(np,:))

         end do
      end do

      j = ny
      do i = 2, nx-1
         np = (j-1)*nx + 1
         nps = np - nx

         cp(np) = cp(nps)
         gcp(np) = gcp(nps)
         Rgp(np) = Rgp(nps)
         hi(np,:) = hi(nps,:)

      end do

      i = 1
      do j = 2, ny-1

         np = (j-1)*nx + i

         ! Information about the P CV to Gibbs module
         call GIBBS_MACH1D_recebe3 (np)

         call GIBBS_EQUILIBRIO_dados_recebe (1, 0, cmod, itimax, itemax, &
            tol_ej, tol_n, T(np), p(np)/101325.0d0, OF)
         if (itt > 1) call GIBBS_MACH1D_recebe6 (ej)

         call GIBBS_EQUILIBRIO_calculos1

         call GIBBS_MACH1D_envia8 (e_j)

         do ii = 1, Nr
            ej (np,ii) = e_j(ii)
         end do

         call GIBBS_hi_massa_calculo
         call GIBBS_MACH1D_envia4 (cp(np), gcp(np), Rgp(np), Yi(np,:), hi(np,:))

      end do

      i = nx
      do j = 2, ny-1

         np = (j-1)*nx + i

         ! Information about the P CV to Gibbs module
         call GIBBS_MACH1D_recebe3 (np)

         call GIBBS_EQUILIBRIO_dados_recebe (1, 0, cmod, itimax, itemax, &
            tol_ej, tol_n, T(np), p(np)/101325.0d0, OF)
         if (itt > 1) call GIBBS_MACH1D_recebe6 (ej)

         call GIBBS_EQUILIBRIO_calculos1

         call GIBBS_MACH1D_envia8 (e_j)

         do ii = 1, Nr
            ej (np,ii) = e_j(ii)
         end do
         call GIBBS_hi_massa_calculo
         call GIBBS_MACH1D_envia4 (cp(np), gcp(np), Rgp(np), Yi(np,:), hi(np,:))

      end do

   end subroutine eq_frozen_cp_gamma

   !****************************************************************************

   ! Subroutine 02

   subroutine eq_effective_gamma(nx, ny, T, p, cmod, itimax, itemax, tol_ej, &
         tol_n, OF, dpY, Rgp, gcpf, gcp) ! The last is output

      integer, intent(in) :: nx     ! Number of control volumes in csi direction (real + ghost)
      integer, intent(in) :: ny     ! Number of control volumes in eta direction (real + ghost)
      integer, intent(in) :: cmod   ! Chemical model
      integer, intent(in) :: itimax ! Maximum number of iterations to solve e(j) - dissociation rate for j-reaction
      real(8), intent(in) :: tol_ej ! Tolerance for e(j) evaluation
      integer, intent(in) :: itemax ! Maximum number of iterations to solve n - the number of mols of products
      real(8), intent(in) :: tol_n  ! Tolerance for n evaluation
      real(8), intent(in) :: OF     ! Oxidant/Fuel ratio
      real(8), intent(in) :: dpY    ! Pressure variation used to estimate gcp

      real(8), dimension(:), intent(in) :: T     ! temperature (K)
      real(8), dimension(:), intent(in) :: p     ! pressure (Pa)
      real(8), dimension(:), intent(in) :: gcpf  ! frozen gamma (non-dimensional)
      real(8), dimension(:), intent(in) :: Rgp   ! gases mixture constant (J/kg.K)
      real(8), dimension(:), intent(inout) :: gcp   ! specific heats ratio (non-dimensional)

      integer :: nps       ! Counter: south
      integer :: i, j, np  ! Counters

      do i = 2, nx-1
         do j = 1, ny

            np = (j-1)*nx + i

            call GIBBS_EQUILIBRIO_dados_recebe (1, 0, cmod, itimax, itemax, &
               tol_ej, tol_n, T(np), (p(np) + dpY*p(np)/2.0d0)/101325.0d0, OF)
            call GIBBS_EQUILIBRIO_calculos1
            call GIBBS_MACH1D_reativo5
            call GIBBS_EQUILIBRIO_dados_recebe (1, 0, cmod, itimax, itemax, &
               tol_ej, tol_n, T(np), (p(np)-dpY*p(np)/2.0d0)/101325.0d0, OF)
            call GIBBS_EQUILIBRIO_calculos1
            call GIBBS_MACH1D_reativo6
            call GIBBS_MACH1D_reativo4 (dpY*p(np), gcpf(np), 8314.51d0/Rgp(np),&
               p(np), gcp(np))

         end do
      end do

      j = ny
      do i = 2, nx-1
         np = (j-1)*nx + 1
         nps = np - nx

         gcp(np) = gcp(nps)

      end do

      i = 1
      do j = 2, ny-1

         np = (j-1)*nx + i

         call GIBBS_EQUILIBRIO_dados_recebe (1, 0, cmod, itimax, itemax, &
            tol_ej, tol_n, T(np), (p(np) + dpY*p(np)/2.0d0)/101325.0d0, OF)
         call GIBBS_EQUILIBRIO_calculos1
         call GIBBS_MACH1D_reativo5
         call GIBBS_EQUILIBRIO_dados_recebe (1, 0, cmod, itimax, itemax, &
            tol_ej, tol_n, T(np), (p(np)-dpY*p(np)/2.0d0)/101325.0d0, OF)
         call GIBBS_EQUILIBRIO_calculos1
         call GIBBS_MACH1D_reativo6
         call GIBBS_MACH1D_reativo4 (dpY*p(np), gcpf(np), 8314.51d0/Rgp(np), &
            p(np), gcp(np))

      end do

      i = nx
      do j = 2, ny-1

         np = (j-1)*nx + i

         call GIBBS_EQUILIBRIO_dados_recebe (1, 0, cmod, itimax, itemax, &
            tol_ej, tol_n, T(np), (p(np) + dpY*p(np)/2.0d0)/101325.0d0, OF)
         call GIBBS_EQUILIBRIO_calculos1
         call GIBBS_MACH1D_reativo5
         call GIBBS_EQUILIBRIO_dados_recebe (1, 0, cmod, itimax, itemax, &
            tol_ej, tol_n, T(np), (p(np)-dpY*p(np)/2.0d0)/101325.0d0, OF)
         call GIBBS_EQUILIBRIO_calculos1
         call GIBBS_MACH1D_reativo6
         call GIBBS_MACH1D_reativo4 (dpY*p(np), gcpf(np), 8314.51d0/Rgp(np), &
            p(np), gcp(np))

      end do

   end subroutine eq_effective_gamma

   !****************************************************************************

   ! Subroutine 03

   subroutine eq_Mach (Rgp, gcp, u, T, M) ! The last one is output
      implicit none
      real(8), dimension(:), intent(in) :: Rgp  ! Gas constant (J/kg.K)
      real(8), dimension(:), intent(in) :: gcp  ! gcp = gamma = Cp/Cv at the center of CV P (non-dim.)
      real(8), dimension(:), intent(in) :: u    ! velocity (m/s)
      real(8), dimension(:), intent(in) :: T    ! temperature (K)
      real(8), dimension(:), intent(out) :: M    ! Mach number

      M = u/dsqrt(gcp*Rgp*T)

   end subroutine eq_Mach

   !****************************************************************************

   ! Subroutine 04

   subroutine eq_frozen_mi_k(nx, ny, T, v_type, visc_cte, k_type, k_cte, vlp, &
         kp) ! The last two are output
      implicit none
      integer, intent(in) :: nx ! Number of volumes in csi direction (real + ghost)
      integer, intent(in) :: ny ! Number of volumes in eta direction (real + ghost)
      integer, intent(in) :: v_type   ! v_type = 1 -> constant; v_type = 2 -> frozen
      real(8), intent(in) :: visc_cte ! constant dynamic viscosity (only for v_type = 1)
      integer, intent(in) :: k_type   ! k_type = 1 -> constant; k_type = 2 -> frozen
      real(8), intent(in) :: k_cte    ! constant thermal conductivity (only for k_type = 1)
      real(8), dimension(:), intent(in) :: T    ! temperature (K)
      real(8), dimension(:), intent(out) :: vlp ! dynamic viscosity (Pa.s)
      real(8), dimension(:), intent(out) :: kp  ! thermal conductivity (W/m.K)

      integer :: i

      if (v_type == 2 .or. k_type == 2) then

         vlp = 0.0d0
         kp = 0.0d0

         do i = 1, nx*ny

            call GIBBS_MACH1D_recebe2 (T(i))

            call GIBBS_CONGELADO_calculos2

            call GIBBS_CONGELADO_envia2 (vlp(i), kp(i))

         end do

      end if

      if (v_type == 1) vlp = visc_cte

      if (k_type == 1) kp = k_cte

   end subroutine eq_frozen_mi_k

   !****************************************************************************

   ! Subroutine 05

   subroutine eq_entropy_and_corners_correction(unit_au, nx, ny, cp_type, &
         ccTw, cmod, Ns, Nr, it, itmax, itimax, tol_ej, itemax, tol_n, po, T0, &
         OF, dpY, Twall, p, pl, M, cp, Rgp, hi, ej, T, gcp, gcpf, Yi, &
         ent, entin) ! Output entries

      implicit none
      integer, intent(in) :: unit_au  ! Auxiliar unit to data storage
      integer, intent(in) :: nx       ! Number of volumes in csi direction (real + ghost)
      integer, intent(in) :: ny       ! Number of volumes in eta direction (real + ghost)
      integer, intent(in) :: cp_type  ! cp_type = 1 -> constant; cp_type > 1 -> variable (evaluated according the temperature/chemical composition)
      integer, intent(in) :: ccTw     ! ccTw = 0 -> adiabatic wall; ccTw = 1 -> prescribed temperature at wall
      integer, intent(in) :: cmod     ! Chemical model
      integer, intent(in) :: Ns       ! Number of chemical species
      integer, intent(in) :: Nr       ! Number of chemical reactions
      integer, intent(in) :: it       ! Number of global iterations
      integer, intent(in) :: itmax    ! Maximum number of global iterations
      integer, intent(in) :: itimax   ! Maximum number of iterations to solve e(j) - dissociation rate for j-reaction
      real(8), intent(in) :: tol_ej   ! Tolerance for e(j) evaluation
      integer, intent(in) :: itemax   ! Maximum number of iterations to solve n - the number of mols of products
      real(8), intent(in) :: tol_n    ! Tolerance for n evaluation
      real(8), intent(in) :: po       ! Stagnation pressure in the chamber (Pa)
      real(8), intent(in) :: T0       ! Stagnation temperature in the chamber (K)
      real(8), intent(in) :: OF       ! Oxidant/Fuel ratio
      real(8), intent(in) :: dpY      ! Pressure variation used to estimate gcp
      real(8), dimension(:), intent(in) :: Twall ! temperature at wall (K)
      real(8), dimension(:), intent(in) :: p     ! pressure (Pa)
      real(8), dimension(:), intent(inout) :: pl    ! pressure correction
      real(8), dimension(:), intent(inout) :: M     ! Mach number
      real(8), dimension(:), intent(inout) :: cp    ! specific heat at constant pressure (J/kg.K)
      real(8), dimension(:), intent(inout) :: Rgp   ! gases mixture constant (J/kg.K)
      real(8), dimension(:,:), intent(inout) :: hi    ! Enthalpy for each chemical species (J/kg.K)
      real(8), dimension(:,:), intent(inout) :: ej    ! Dissociation rate
      real(8), dimension(:), intent(inout) :: T     ! temperature (K)
      real(8), dimension(:), intent(inout) :: gcp   ! specific heats ratio (non-dimensional)
      real(8), dimension(:), intent(inout) :: gcpf  ! frozen gamma (non-dimensional)
      real(8), dimension(:,:), intent(inout) :: Yi    ! Mass-fraction of chemical species (non-dimensional)
      real(8), dimension(:), intent(inout) :: ent   ! Entropy in each control volume(J/kg.K)
      real(8), dimension(:), intent(inout) :: entin ! Entropy at the chamber inlet (J/kg.K)

      integer :: i, j, ii  ! Counters
      integer :: np, nps, npn, npe, npw   ! Counters - control volumes

      ! Storage of auxiliary fields
      if (cp_type > 1) then
         ! Auxiliar resource to avoid negative temperatures at wall (Part 01)
         if (ccTw == 1) then
            j = ny
            do i = 1, nx
               np = (j-1)*nx + i
               T(np) = Twall(i)
            end do
         end if
         call GIBBS_EQUILIBRIO_dados_recebe (unit_au, 0, cmod, itimax, itemax, &
            tol_ej, tol_n, T0, po/101325.0d0, OF)
         call GIBBS_EQUILIBRIO_calculos1
         call GIBBS_si_aloca_memoria
         call GIBBS_MACH1D_camara_entropia (entin, ny)
         call GIBBS_EQUILIBRIO_dados_escreve (unit_au)
         call GIBBS_EQUILIBRIO_escreve1
         call GIBBS_CONGELADO_calculos2
         call GIBBS_CONGELADO_escreve2 (unit_au)
         call eq_frozen_cp_gamma(nx, ny, T, p, cmod, Nr, it, itmax, itimax, &
            itemax, tol_ej, tol_n, OF, Yi, hi, Rgp, cp, gcp, ej)
         gcpf = gcp
         call eq_effective_gamma(nx, ny, T, p, cmod, itimax, itemax, tol_ej, &
            tol_n, OF, dpY, Rgp, gcpf, gcp)

         do i = 2, nx-1
            do j = 1, ny
               np = (j-1)*nx + i
               call GIBBS_MACH1D_entropia (T(np), p(np)/101325.0d0, Rgp(np), &
                  Yi(np,:), ent(np))
            end do
         end do

         i = 1
         do j = 2, ny-1
            np = (j-1)*nx + i
            call GIBBS_MACH1D_entropia (T(np), p(np)/101325.0d0, Rgp(np), &
               Yi(np,:), ent(np))
         end do

         i = nx
         do j = 2, ny-1
            np = (j-1)*nx + i
            call GIBBS_MACH1D_entropia (T(np), p(np)/101325.0d0, Rgp(np), &
               Yi(np,:), ent(np))
         end do

      end if

      ! Corrections at corners
      ! Corner: SW (lower-left)
      j = 1
      i = 1
      np = (j-1)*nx + i
      npn = np + nx
      M(np) = M(npn)
      gcp(np) = gcp(npn)
      Rgp(np) = Rgp(npn)
      ent(np) = ent(npn)
      pl(np) = 0.0d0
      do ii = 1, Ns
         Yi(np,ii) = Yi(npn,ii)
      end do

      ! Corner: SE (lower-right)
      j = 1
      i = nx
      np = (j-1)*nx + i
      npn = np + nx
      M(np) = M(npn)
      gcp(np) = gcp(npn)
      Rgp(np) = Rgp(npn)
      ent(np) = ent(npn)
      pl(np) = 0.0d0
      do ii = 1, Ns
         Yi(np,ii) = Yi(npn,ii)
      end do

      ! Corner: NW (upper-left)
      j = ny
      i = 1
      np = (j-1)*nx + i
      nps = np - nx
      npe = np + 1
      M(np) = M(npe)
      gcp(np) = gcp(npe)
      Rgp(np) = Rgp(npe)
      ent(np) = ent(npe)
      pl(np) = 0.0d0
      do ii = 1, Ns
         if (ccTw == 0) then
            Yi(np,ii) = Yi(nps,ii)
         else
            Yi(np,ii) = Yi(npe,ii)
         end if
      end do

      ! Corner: NE (upper-right)
      j = ny
      i = nx
      np = (j-1)*nx + i
      nps = np - nx
      npw = np - 1
      M(np) = M(npw)
      gcp(np) = gcp(npw)
      Rgp(np) = Rgp(npw)
      ent(np) = ent(npw)
      pl(np) = 0.0d0
      do ii = 1, Ns
         if (ccTw == 0) then
            Yi(np,ii) = Yi(nps,ii)
         else
            Yi(np,ii) = Yi(npw,ii)
         end if
      end do

   end subroutine eq_entropy_and_corners_correction

   !****************************************************************************

end module equilibrium
